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Richardson equations can be mapped on the classical electrostatic problem in two di- 
mensions. We have recently suggested a new analytical approach to these equations in the 
thermodynamical limit, which is based on the 'probability' of the system of charges to be in 
a given configuration at the effective temperature equal to the interaction constant. In the 
present paper, we apply this approach to excited states of the Richardson pairing model. 
We focus on the equally-spaced situation and address arbitrary fillings of the energy layer, 
where interaction acts. The 'partition function' for the classical problem on the plane, which 
is given by Selberg-type integral, is evaluated exactly. Three regimes for the energy gap are 
identified, which can be treated as the dilute regime of pairs, BCS regime, and dilute regime 
of holes. 

Subject Index: 010, 062, 368 

§1. Introduction 

Bardeen-Cooper-Schrieffer (BCS) theory plays a very important role in the mi- 
croscopic description of superconductivity.®'®'® As shown by Richardson,® BCS 
Hamiltonian turns out to be exactly solvable. By staying in the canonical ensemble, 
Richardson managed to find a many-body wave function, which depends on the set of 
energy-like quantities (rapidities) and provides an exact solution of the Schrodinger 
equation. The number of rapidities is equal to the number of Cooper pairs in the 
system, while the Hamiltonian eigenvalue is given by their sum. Rapidities satisfy 
the system of nonlinear algebraic equations, now called Richardson equations. The 
resolution of these equations is a formidable task. More recently, it was shown that 
Richardson equations can be derived from the algebraic Bethe-ansatz approach.® 
They are also closely related to the well-known Gaudin mode® and Chern-Simons 
theories?' Richardson equations are now widely used to study numerically super- 
conducting state in nanometer-scale systems.® 1 ® In particular, powerful tools of the 
quantum inverse scattering method were used to compute correlation functions^® in 
such systems. 

It is quite remarkable that Richardson equations can be mapped onto the clas- 
sical electrostatic problem in the plane.®'^ Namely, energy-like quantities may 
be treated as coordinates of interacting charged particles, which are placed into the 
external electric field. Equilibrium positions of these charges are then equivalent to 
solutions of Richardson equations. The origin of this highly remarkable example of 
quantum-to-classical correspondence is unclear. 

We recently pushed further^® the analogy with the classical electrostatic prob- 
lem by introducing the occupation 'probability' for the system of charges at the 
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effective 'temperature' given by the interaction amplitude, which goes to zero in 
the thermodynamical limit. This leads to the new approach to treat analytically 
Richardson equations. Namely, one can reconstruct an information on the location 
of the center of masses for the free charges by using an integration instead of a 
straightforward resolving the equations. This is done by constructing the 'partition 
function', given by the Selberg-type integral, so that the energy of the initial quantum 
problem is determined by the logarithmic derivative of this classical quantity with 
respect to the inverse temperature. Note that Selberg integrals are familiar in con- 
formal field theory and in random-matrix models. This fact suggests an interesting 
link with these subjects. 

In the present paper, we extend the 'probabilistic' approach to excited states 
of Richardson pairing model, which amounts manipulating more sophisticated and 
general Selberg-type integrals. The 'partition function' for the deterministic problem 
is found analytically by converting such an integral into a coupled binomial sum, 
evaluated using combinatorial properties of Vandermonde matrix. We then calculate 
energy difference (gap) between excited and ground states, which is a more subtle 
quantity than the ground state energy itself. At the same time, this quantity is much 
more important, since it is not possible to test ground state energy experimentally 
in contrast to the gap. We focus on the equally-spaced model and treat arbitrary 
fillings of the energy interval, where attraction between up and down spin electrons 
acts. This can be seenPJ as a toy model for the density-induced crossover between 
individual fermionic molecules and a dense regime of BCS pairsP^S We identify 
three different regimes for the energy of the excited state. The first one corresponds 
to low densities of pairs; the excitation energy in this case is controlled by the single- 
pair binding energy. The second regime corresponds to the dense regime of pairs. 
It is described by the gap of BCS type, which has a cooperative origin. The third 
regime is a 'superdense' regime of pairs or a dilute regime of holes, which is again 
controlled by the single-pair binding energy. 

The paper is organized as follows. In §2, we briefly formulate the problem and 
outline basic ingredients of our method. In §3, we calculate the 'partition function'. 
In §4, we discuss our results and we conclude in §5. 

§2. General formulation 

We consider electrons with up and down spins, which interact through the BCS 
potential. The Hamiltonian is given by 

H = Yl £k (4t a kt + a U a n) ~ V Y1 4't a -k'4 a -H«kt- (2-i) 

k k,k' 

We stay in the canonical ensemble, i.e., with the number of electrons fixed. In 
this case, the eigenvalue of Hamiltonian, given by equation (|2-ip . can be represented 
as a sum of N energy-like quantities Rj, which satisfy the system of ./V coupled 
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Richardson equations 

It is assumed that the attraction between spin up and spin down electrons acts 
only for the free electron states having kinetic energies confined between ep and 
£f + 17; the same therefore applies to the sums in the right-hand side (RHS) of 
equation (|2-2|) . The first quantity ep , within the traditional BCS framework, can be 
associated with the Fermi energy of the filled Fermi sea of noninteracting electrons, 
while 17/2 is the Debye frequency. We then also assume that all the energies of 
free electron states are different and they are distributed equidistantly within this 
potential layer. Such a situation is usually refereed to as the equally-spaced model. 
It can be understood in terms of a constant density of energy states. The distance 
between two nearest levels is given by 1/p, where p is the density of energy states. 
This means that we consider a situation with the discrete energy spectrum, but the 
distance between nearest energy levels for noninteracting electrons goes to zero in 
the thermodynamical limit, so that it becomes quasi-continuous. Thus, there are in 
total Nq = pQ available energy states of each spin direction in the potential layer. 
We fill this potential layer by N pairs. In the usual BCS configuration, N = Nq/2 
(half-filling) , while we address the arbitrary filling corresponding to macroscopic N 
(finite density of pairs). In other words, we consider a thermodynamical limit, for 
which the dimensionless interaction constant v = pV stays independent on N, as 
well as the filling N/Nq, whereas p ~ N . Now we take the limit N — > oo. Changing 
the filling should be considered as a toy model for the density-induced crossover from 
the dilute regime of pairs, when the filling is small, to the dense regime, when the 
filling is increased PIiJ<G3B This crossover attracts a lot of attention in the field of 
ultracold gases. In addition, it can be relevant for high-T c cupratesP^'^ Such a toy 
model also helps to establish a link between the single-pair problem solved by Cooper 
and many-pair BCS condensate. It was arguecP^ that this model can be directly 
applicable to some semiconductors. Moreover, it has obvious similarities with the 
well-known Eagles model,--' which however does not assume constant density of 
states. 

There are two types of lowest excited states in Richardson approach to the pair- 
ing problem. Excitations of the first kind correspond to one of the rapidities having 
zero imaginary part and located between two free electron levels, i.e., in a quasi- 
continuum spectrum. Excitations of the second type correspond to the presence 
of unpaired electrons, which can appear for instance due to the breaking of one of 
the pairs. These unpaired electrons then block the one-electronic states, which they 
populated Their role is to bring their own bare kinetic energy to the total energy of 
the system, and, moreover, to modify the energy of the remaining set of pairs. Here 
we consider this second type of excited states by addressing a system of N pairs plus 
one unpaired electron. At the end, we will show how, in the thermodynamical limit, 
the energy of the excited state of the first kind can be found by reducing the problem 
to the excited state of the second kind. Note that, within the variational approach 
of Ref. [T]), excitations of the first type were called 'real pairs', while excitations 
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of the second type correspond to 'broken pairs'. Within the Bogoliubov approach 
both types of excitations are handled on the same footing,® without making any 
distinction between them. 

The role of a single unpaired electron is to block the occupied state: No scattering 
occurs, since BCS interaction potential couples only electrons with opposite spins 
and momenta.®'® Hence, this state must be excluded from the sums appearing in 
Richardson equations (|2-2p . We denote the total lowest possible energy of the system 
of N pairs plus one unpaired electron with momentum p as E^ p . For the energy 
of N pairs in presence of the unpaired electron we use Er p . The latter quantity 
can be found from Richardson equations, where the corresponding blocked state is 
excluded from the sums. Consequently, En )P = -Er, p + £ p - The ground state energy 
of N pairs is denoted as Ejf. Within the 'probabilistic' approach, En was calculated 
in Ref. Q2D. 

It was shown long time ago that Richardson equations have a remarkable elec- 
trostatic analogy.^ '® Consider the function E c [ ass ({Rj}), given by 

E c iass{{Rj}) = 2 \ Y1 ReR i + V Y ln l 2ek - R i\- 2V Y ln \ R i - R A> ( 2 ' 3 ) 

\ 3 j,k 3,l,j<l j 

which can be interpreted as the energy of ./V free classical particles with electrical 
charges 2y/V located on the plane with coordinates given by (Re Rj, Im Rj). Free 
particles are placed into a uniform external electric field, which gives rise to the force 
(—2, 0) acting on each particle. In addition, free charges repeal each other, and they 
are also attracted to Nq fixed particles each having a charge — y/V and located at 
(2£k)0). Richardson equations can be formally written as the equilibrium condition 
for the system of N free charges. This can be seen by splitting E c i ass ({Rj}) as 

W({Rj}) + W(|i?||) with 

W({Rj}) =J2 R j + v J2 ln ( 2ek " R ^ ~ 2V Yl ln - R ^ ( 24 ) 

j i,k j,i,j<i 

and by considering conditions dW({Rj})/dRj = 0, from which Richardson equations 
(|2-2p follow. Note that the equilibrium of this system of free charges is not stable. If 
we treat the excited state with one level blocked, the corresponding term must be, 
of course, excluded from the sums in the RHS of equations (|2-3p and (|2-4p . while, in 
the case of the ground state, the summation runs over all terms. 

We recently suggested an idea to push further the analogy between the initial 
quantum problem and the classical problem of Coulomb plasma in two dimensions 
by considering occupation 'probabilities' S({Rj}) = exp (— W({Rj} /T e ff) at the ef- 
fective 'temperature' T e ff = V. The effective 'temperature', in the thermodynamical 
limit, goes to zero as 1/N, so one can reconstruct an information about the sum of 
energy-like quantities in equilibrium without solving Richardson equations directly, 
but using an integration of S({Rj}) over Rj, in a spirit of usual thermodynamics, 
except of the facts that integration is performed over half of degrees of freedom 
(one-dimensional integration over each Rj) and S({Rj}) is a meromorphic function, 
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which allows one to deform integration paths. The approach has some similarities 
with the large- N expansion for the Dyson gasA^* 

In the case of one level blocked, taking into account an above expression of 
S({Rj}), we write it in the form 



S({Rj}) 



exp 



(2-5) 



where the blocked level 2e p = 2ef + 2n p /p has been excluded from the product in 
the denominator. The 'probability' can be represented as 



S({Rj}) 



niiiikC^k-fl,-) 



exp 



V 



(2-6) 



Note that, while S({Rj}) for the ground state has analogies with the square of 
Laughlin wave function for the ground statep^ equation (|2-6|) resembles Laughlin 
wave function for excited states, since it has an additional factor of the same type, 
nj(2e p — Rj), compared to the ground-state S. 

Next, we perform a partial- fraction decomposition of S and rewrite it as 



S({Rj = 2e Fo - rj}) = exp 



-2Ne Fo 
V 




if , 2n p- 



II (n - r 3 f 

},3,l>3 
N n \ exp(r i /y) 



n ; 



r,- H - 

3 p 



(2-7) 



where is the binomial coefficient. Here and below we drop all irrelevant pref- 

\ rij / 

actors independent on V. 

To apply our technique, we must be sure that S({Rj}) goes to zero, as an 
imaginary part of any of the R tends to infinity. This criterium is satisfied provided 
that the filling is not larger than 1/2, as follows from equation f|2-6|) . To address 
larger fillings, we should switch to the representation in terms of holes, as it has been 
done in Ref. fT2l). 



§3. Evaluation of 'partition function' 



3.1. 'Partition function' as the binomial sum 

We now introduce a 'partition function', which is given by the Selberg-type 
multidimensional integral 

Z = j S({Rj})dRj, (3-1) 

where an integration is performed over the whole set {Rj} along the contour for each 
Rj, shown in Fig. 1. The integration path avoids all the poles of S (corresponding to 
2ffk) and then reconnects via the semicircle of infinite radius, along which S({Rj}) is 
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zero. The sum -E_r, p of energy-like quantities R in the equilibrium, which is a quantity 
of interest, is expressed through the logarithmic derivative of Z with respect to the 
inverse 'temperature' V 

E R>p = --^hxZ. (3-2) 




Fig. 1. The schematic plot of complex values of each Rj. Filled circles show locations of energy 
levels for free electrons. Open circle indicates a missing (blocked) level. Solid line corresponds 
to the integration path for each Rj. 

At this stage, we substitute equation (|2-7p to equation (|3-1|) and integrate using 
residues. We easily get 

Z = exp(-2Ne Fo /V)z, (3-3) 

where 




where a = exp(—2/v). 
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3.2. Vandermonde determinants 



It is known that Y\.iji>j{ n l ~~ n j) entering equation ([3 
the determinant of the Vandermonde matrix 



Y[ ( n i - n j) = v o({ n j}) = det 

i,j,i>j 



n 



1 


1 


1 


1 


ni 


n 2 


n 3 




n\ 


nl 


nj 


• < 


N-l 
1 




nf 1 : 





can be represented as 



(3-5) 



Consequently, we can rewrite the total product Y\i j i>j(. n i ~ n j) 2 as 
(-l) N ( N - 1)/2 V ({ nj })Vo{{N a - nj }). 

Instead of the standard monomials, we represent Vo({nj}) in terms of Pochham- 
mer symbols (n) a = n(n — l)...(n — a + 1), while (n)o = 1. This can be done by 
iterative substraction of rows of Vandermonde matrix, which does not change the 
determinant and ultimately leads to 



V ({ nj }) = det 



( n i)i 



(n 2 )o 
(n 2 )i 
(n 2 ) 2 



Ml 

(n 3 ) 2 



(ni) 



N-l 



(«2) 



N-l 



(ns) 



N-l 



(n N ) 

(njv)i 
(n N ) 2 



(nN) 



N-l 



(3-6) 



while Vq({Nq — n>j}) can be represented in a similar form with nj changed into 
Nq - nj. 

Next, we insert (nj—n p ), which appears in equation (|3-4p . into the Vandermonde 
matrix in the following manner 



N 



det 



(ni)o(ni - n p ) 
(ni)i(ni - rip) 



(n 2 ) (n 2 - rip) 
(n 2 )i(n 2 - n p ) 




(ni)jv-i(ni -rip) (n 2 )Ar-i(n 2 - n p ) 



n ("i - n J 

3,l,l>3 

(n N ) (n N - rip) 
(nAr)i(njv - n p ) 

(nAr)Ar_i(nAr - n p ) 



(3-7) 



Let us now extract n p from the first row of the above matrix as 







("2)1 




(njv)i 




det 


(ni)i(ni - n p ) 


{n 2 )i(n 2 - 


Hp) 


(nAr)i(nAr- 


-P) 




.(«l)iV-i(ni - n p ) 


(n 2 )N-i(n 2 


-Up) . 


• • (nN)N-i(nN 


-n p )_ 




(«i)o 


(n 2 )o 




("3)0 




n p det 


(ni)i(ni - n p ) 


(n 2 )i(n 2 - 


n p ) 


(nAr)i(nAT- 


n p) 




_(ni)Ar„i(ni - n p ) 


(n 2 ) N -i(n 2 


-Up) . 


•• (n N ) N -i(riN 


-n p )_ 



= £>1 



n p D 2 . 



(3-8) 



8 



W. V. Pogosov 



The first term, D\, can be rewritten as 



det 



1 det 



(«l)l 
("1)2 

(ni)jv-i(ni - rip) 

( n i)i 
(«i)l 



("2)1 
("-2)2 

(n 2 )jv-i(n 2 - n p ) 

("•2)1 
(n 2 )i 



(ni)Ar-i(ni - n p ) (n 2 )jv-i(n 2 - n p ) 



[113)1 
(n N ) 2 

(n N ) N _i(n N - n p ) 

(njv)i 
(njv)i 



(n N ) N -i(n N - n 



(3-9) 



where we have extracted ra p — 1 from the second row of the initial matrix. We im- 
mediately see that the last determinant is exactly zero, since two rows of its matrix 
are the same. We then extract n p — 2 from the third row of the remaining ma- 
trix in the RHS of equation (|3-10p and repeat our arguments. We follow a similar 
iterative procedure for D 2 . It is rather straightforward to see that only terms pro- 
portional to (n p ) m do survive at the end, since all the other terms are proportional 
to determinants of matrices with repeating rows. We finally arrive to the following 
identity 



N 

n< 



11; 




N 



^(-l) m K) m %- m (K}) (3-10) 



m=0 



where V m ({n,j}) are determinants of the matrices, which are obtained from the matrix 
of the RHS of equation (|3-6p by increasing indices of Pochhammer symbols in the 
last m rows by one, as (rij) a — > {nj) a+ \. 

3.3. Auxiliary identities 

To go further, we provide some auxiliary identities, which are going to greatly 
facilitate calculations. The starting identity is 



Za,b 



N„ 

E(-i: 

n=0 



n a n ( Nn )(n) a (Na-n) b 
n 



o a (l - a 



.Nn-a-b 



(-1)° 



(N, 



n 



a -by: 

(3-11) 

where a + b < Nq. It can be obtained by observation that first a terms, as well as 
last b terms of the initial sum are zero, and then by replacing Pochhammer symbols 
by ratios of factorials. 

Let us now focus on a product of sums for each rij, every sum being similar to 
the one of equation ()3Tip . We can write such a product as 



N 

n 



Zai ,bi • • "^ajv ,6jv 

[{Nn-aj-bjY' 



(3-12) 
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dependencies on sets a,j and bj is only by ^2f=\ cij and Ylf=i bj, which actually are 
degrees of polynomials Ylf =1 {nj) aj and Ylf =l (Na - rij) bj , respectively. 



The dependence of this quantity on V (via a) is through two first factors. Their 

Jj=1 dj and J2 j= 

3.4. 'Partition function': hyper geometric series 

Let us substitute equation (|3-10p back to equation (|3-4p . To proceed in calcula- 
tions, we have to consider products of the form V m ({rij})Vo({Nn — nj}). It is easy 
to see from equation (|3-6p that each of them (for a given m) can be represented as 

a linear combination of polynomials of the form (j^\jLi(nj) aj S j (Y\f=i(^^ ~ n j)b^j 

with the unique Y$=i a j = m + N ( N ~ l )I 2 and Ejli b j = N(N - l)/2 for each 
polynomial. These two numbers just give the degrees of the polynomials V m ({nj}) 
and Vo({Nfi — nj}), respectively. According to equation (|3-12p . after substitution 
of these products to equation (|3-4p and performing summations, we get the same 
dependence of the result on a for any of these polynomials (at a given m). Hence, 
we can write 

N 

z = a N(N + l)/2 (1 _ a) N(N n -N) ^ K)m (o ._i _ ^ 

m=0 

where a m are unknown numbers independent on a, which have a combinatorial 
origin. 

We avoid a direct calculation of a m by using a trick, which is based on the well- 
known rule: ( N °) = (^J . Namely, we change summation variables in equation 

(|3-4p as rij = Nq — rij. It is then readily seen that 

z(n p ,a) = (-l^+^a^ziNn-n^a- 1 ). (3-14) 

For z{Nq — ?ip, cr _1 ) we can use equation (|3T3p with n p — > — n p , a — > a^ 1 . By 
doing this, we arrive to another expression of z(n p , a) in terms of (Nq 

— ^p)m; which 

is nevertheless equivalent to equation (|3-13p : 



TV 

z = a N(N-i)/2 (1 _ a) N( Nn -N) {Nn - n p ) m (a" 1 - l) m a m . (3-15) 

m=0 

The idea is to express Pochhammer symbols of Nq — n p in terms of Pochhammer 
symbols of n p , then to substitute them to equation (|3-15p and to compare the result 
with equation (|3T3p . By equating coefficients of Pochhammer symbols of n p , we are 
going to obtain a system of linear equations for {a m }- We make use of the following 
relation 

(^-n^B-!)^)^^, (3,6) 

which can be trivially checked for m = 0, 1 and then proved by induction. Inserting 
it to equation (|3-15p and solving the system of equations for {a m }, we get 

a m = a() \ , (3-17) 
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where a is an irrelevant constant, independent on a. 
Finally, we obtain z as 



z = aa N ( N+i y 2 (l-af^~ N h', 



(3-18) 



where 



z 




ni 



(3-19) 



The last expression can be considered as a hypergeometric series. 

3.5. Hypergeometric series: saddle-point method 

Actually, equation <|3- 18|> already allows us to find Er p in terms of the hyper- 
geometric series. The resulting expression, however, is essentially untractable from 
the perspective of a further analysis. Let us, therefore, try to transform z' into a 
simpler form. 

We note that, for natural numbers n p and Tjij {jip)m — for any uij which 
is larger than n p . Hence, we can change the upper limit of summation in equation 
(|3-19j> to min(n p , N). After that, all the terms in the sum are nonzero, so that we can 
replace Pochhammer symbols of n p by ratios of factorials, as (n p ) m = n p !/(n p — m)\. 

We also see that a -1 — 1 = exp(2/w) — 1 is always positive. The last circum- 
stance is very important: it means that terms in the sum are not oscillating in sign. 
Therefore, we can switch from summation to integration. We also utilize asymptotic 
expansion for factorials entering both Pochhammer symbols and the binomial coeffi- 
cient, since we are interested in macroscopic numbers, ~ N. Note that the important 
case of n p = must be considered as n p = x N, although this particular situation 
can be analyzed separately without switching to integration, since the derivation 
turns out to be quite simple (the results of both approaches are finally the same). 
After straightforward algebra, we arrive to the identity 



To avoid possible confusion, we note that in equation (|3-20p we have dropped, as 
usual, an irrelevant ^-independent prefactor. 

In order to evaluate the integral in the RHS of equation (|3-2U|) . we use a saddle- 
point method. It is straightforward to prove that H{m) has a single maximum within 
the integration range, which is attained at 




(3-20) 



where 



H(m) = {Nq — m) \b.{Nq — m) — mln(m) — [N — m) ln(iV — m) 
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We can also ensure that the position of this maximum is far enough from both 
integration limits, since H"(mo) ~ A -1 , while l/\/ H"(mo) ~ y/~N determines the 
width of the neighborhood of mo, which gives the dominant contribution to the 
integral. At the same time, both mo and (min(n p , A) — tuq) do scale as A. Together 
with the fact that H^ n \mo) ~ N~ n+1 at n > 2, this enables us to reduce the problem 
to the simple Gaussian integration. Keeping leading order in A, we obtain 

lnz' = H(m ). (3-23) 

By finding a logarithmic derivative of z, given by equation (|3-18j> . and by adding 
the kinetic energy of the unpaired electron spo + n p /p, we arrive to the expression 
of the total energy -E/v,p of the system of A pairs and one unpaired electron. After 
some simple algebra and using the known expression for the ground state energy^ 
En, we can present E^, P as 

E N , P = E N + p + ^A 2 + (e p -p) 2 , (3-24) 

while 

A = --^-^N(N n -N), (3-25) 
p 1 — a 

A 1 + a „ a , 

P = e F0 + -- n- . (3-26) 

p 1 — a 1 — a 

Physical meanings of A and p will be fixed below. 

The formalism presented above is restricted to < Nq/2. In order to address 

configurations with N > Nq/2, we switch to holes. It is then straightforward to 
ensure that equation (|3-24p still holds in this case. 

§4. Discussion 

We now consider Ejy iP as a function of e p . In particular, it is of interest to 
determine the lowest possible -EW iP . 

First of all, we see that p, given by equation (|3-26|) . can be actually considered 
as a chemical potential. Indeed, by calculating (En — En-i)/2, we find that this 
quantity does coincide with p. 

The minimum value of -Ejv.p is attained at n p , which gives a minimum of (e p — 
p) 2 , provided that n p is confined between and Nq. If p — sfo also falls into this 
range, one can always choose n p /p = p — efo, so that (e p — p) 2 is zero, while the 
square root in the RHS of equation ()3-24p reduces to A. It is easy to see that for 
the half-filling configuration, Nq = N/2, the expression of A, as given by equation 
(|3-25p . reproduces precisely the BCS formula for the gap. This conclusion is also in 
agreement with the result of the traditional method to solve Richardson equations 
in the large-sample limit, which uses an assumption that energy-like quantities are 
arranged into arcs on the complex plane JIIJ'GSliJZIJ The assumption on arcs is actually 
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deduced from numerical solutions of Richardson equations, so that, in general case, 
it is not fully controllable, from our point of view. Similar approach is also widely 
used in a broader context for the solution of the Bethe-ansatz equations Another 
method,^ 1 which is based on Taylor expansions of sums appearing in Richardson 
equations around the known single-pair solution, up to now has been successfully 
applied to the ground state only, while its application to excited states leads to 
heavy mathematical problems. 

If the chemical potential goes below the lower cutoff epo, i.e., (jl — efo becomes 
negative, a constrained minimum of E^,p corresponds to n p = 0. By solving the 
equation /i — epo = 0, we see that the transition between the two regimes occurs at 
N = No, where 

N = Naj^. (4-1) 

If the chemical potential goes above the upper limit of the potential layer, the 
minimum energy corresponds to e p located at this upper limit. By solving the 
equation {jl — efo = ^> we find that the transition to this regime happens at Nq — Nq. 
This result is in agreement with the electron-hole symmetry. 

Thus, for the minimum of E^,p, we have three regimes 

min(£jv iP ) =E N + n+Hj + ' ^ 

at N < N ; 

mm(E N>p ) = E N + fi + A, (4-3) 



at N < N < N n - N 



mm(E N)P ) = E N + p + j (il - + ^3-3^} ' 



(4-4) 



at N n -N < N. 

In the weak-coupling limit, a <C 1, we have: iVo ~ Nqo. This value actually 
corresponds to the density of pairs, at which their wave functions start to overlap. 
Such an overlap can be considered as a signature of the transition from the isolated- 
pair regime to the dense condensate. 

In the dilute regime, iV < iVo, the excitation energy is controlled by Qa / (1 — a), 
which is nothing but half the binding energy e c of an isolated pair. 13 ' In the dense 
regime, iV^ — iVo > iV > iVo, it is governed by the BCS gap A, which has to be 
considered as a collective many-body response of the system to the appearance of 
one blocked level. It is of interest to note that pair binding energy and the gap have 
similar, but different dependencies on interaction constant v. Namely, in the weak- 
coupling limit, the first quantity is proportional to exp(—2/v), while the second one 
behaves as ~ exp(—l/v) due to equation (|3-25p . This equation also shows that A 
is symmetric with respect to the mutual replacement of electron and holes, so that 
the electron-hole symmetry again shows up. The third regime can be considered 



Excited states in Richardson model 



13 



as a 'superdense' regime of Cooper pairs, made of electrons, or the dilute regime 
of Cooper pairs, made of holes. In this regime, we again see a single-pair binding 
energy appearing in min(£7v p ). It can be now understood as a binding energy of a 
pair made out of holes. 

Thus, we have identified three regimes for the energy difference between Ew tP 
and En depending on the energy layer filling. These are a dilute regime of pairs, BCS 
regime, and a dilute regime of holes. With changing layer filling, transitions between 
these regimes occur smoothly. We have found that only higher-order derivatives of 
En,p with respect to N experience discontinuities upon the transitions. 

Up to now, we considered states with only one unpaired electron. Let us address 
the energy of the state with two such electrons. This should allow us to make a com- 
parison with the energy of the system with all the electrons paired, the total number 
of particles being conserved. In principle, in order to find the energy, we should 
perform similar calculations, but with two states blocked. However, we may use a 
simple trick allowing one to avoid making such computations. It is rather obvious 
that it is energetically favorable for the two unpaired electrons to occupy two neigh- 
boring states rather than to be separated. We then come back to the electrostatic 
picture and perform a coarse-graining of the initial configuration. Namely, we merge 
couples of neighboring fixed particles, as well as couples of neighboring free particles 
into 'superparticles' with charges being twice larger than initial ones. We must also 
increase by a factor of two a homogeneous forces acting on free charges. This proce- 
dure does not change dominant (extensive) contribution of the total energy. Within 
this procedure, two blocked states are converted into a single one, so that we can 
map this picture to the above situation. By doing this, after some straightforward 
calculations, we finally reach an expectable result: the difference of energies of the 
states with N pairs and N — 1 pairs plus two unpaired electrons is twice the square 
root entering RHS of (EF22j) . 

Finally, let us discuss another type of excited states, which corresponds to one 
of the energy- like quantities trapped between two one-electronic levels. The energy 
of such a state can be again handled using electrostatic analogy without making 
detailed calculations. Namely, we note that the role of the free charge trapped is 
to compensate two fixed charges, between which it is accommodated, since each of 
them has an opposite charge, but twice smaller in absolute value. We again map this 
configuration to the previous one with two states blocked and find the same gap, 
within dominant terms in 1/N. 

§5. Conclusions 

Richardson equations provide an exact solution for the BCS pairing Hamilto- 
nian. These equations are deterministic and posses a well-known electrostatic anal- 
ogy. Therefore, one can convert the problem of their resolution to the probabilistic 
problem, as it was recently suggested in Ref. 12) for the ground state of the initial 
quantum problem. This approach avoids assumption that Richardson solutions in 
the large- ./V limit are arranged in arcs on the complex plane. 

In the present paper, we applied this treatment to excited states with the fo- 
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cus on the equally-spaced model and the thermodynamical limit. We have con- 
sidered arbitrary fillings of the energy interval, where the attractive potential acts. 
The 'partition function' for the deterministic problem has been found analytically 
by converting the Selberg-type integral into coupled binomial sum, evaluated using 
combinatorial properties of Vandermonde matrix. 

For the energy difference between the first excited state and the ground state 
(energy gap), three regimes have been identified, which can be considered as the 
dilute regime of pairs, BCS regime, and dilute regime of holes. Explicit expressions 
have been derived. Transitions between these regimes occur smoothly, accompanied 
by only weak singularities. The results supports the BCS result for the half-filling. 
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